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1. Introduction 

Why should we use lattice Monte Carlo methods to describe a real sys- 
tem when there are so many efficient and accurate methods to treat the 
full problem? There are at least two good reasons. The first one takes a 
pragmatic point of view: For complicated systems a full calculation simply 
cannot be done on present-day computers. The second reason rests on the 
belief that the physics underlying the properties of real materials is simple 
and can be captured in model systems. If we succeed in this there is the 
additional benefit of having identified the important general features of the 
material. 

For the doped Fullerides we encounter just such a situation. Even for 
a single Cqq molecule a full QMC calculation is still a challenge, and sim- 
ulations of Fullerides, i.e. solids made of Ceo molecules, are simply out 
of question. For many properties it is, however, sufficient to focus on the 
valence band only, removing all other degrees of freedom from the Hamilto- 
nian. Important features of the doped Fullerides that have to be reflected in 
such a model are the degeneracy of the molecular orbital that gives rise to 
the valence band, the filling of the valence band, and the lattice structure 
of the solid. All these can be incorporated in a Hubbard-like Hamiltonian, 
which can be treated efficiently using lattice QMC methods. 

In the following we will first show how to set up a model Hamiltonian for 
the doped Fullerides. Then we discuss Monte Carlo methods for such lat- 
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tice Hamiltonians, especially the optimization of Gutzwiller functions both 
in variational and fixed-node diffusion Monte Carlo. Finally we use QMC 
to investigate the Mott transition in the doped Fullerides. The interest in 
these questions comes from the following situation: Density functional cal- 
culations predict that the doped Fullerides are metals. On the other hand, 
one finds that the Coulomb repulsion between electrons on the same Ceo 
molecule is very strong. This suggests that correlations should be dominat- 
ing, making all doped Fullerides Mott insulators. Reality falls in between 
these two extremes: some doped Fullerides are metals (and even supercon- 
ductors), while others are insulators. From our QMC calculations we find 
that due to the degeneracy of the valence band the integer-doped Fullerides 
are close to a Mott transition, and not far into the Mott insulator regime, 
as simple theories would suggest. Whether a given compound is on the 
metallic or the insulating side of the transition depends then on the crystal 
structure (bipartite vs. frustrated) and the filling of the band. 



2. Model Hamiltonian 

Solid Ceo is characterized by a very weak 
inter-molecular interaction. Therefore the discrete 
molecular levels merely broaden into narrow, well 
separated bands (see Fig. 1) [1]. The valence band 

originates from the lowest unoccupied molecular 

orbital, which is a 3-fold degenerate t\ u orbital. 
Doping the solid with alkali metals does not affect 
the band structure close to the Fermi level very 
much. Only the filling of the t\ u band changes, 
since each alkali atom donates its valence electron. 

To simplify the description of the doped Fullerides 

we want to focus on the electrons in the t\ u band |-| 

only. To get rid of the other degrees of freedom u 

we use the Lowdin downfolding technique [2]. The Figure 1. ■. Schematic band 
basic idea is to partition the Hilbert space into a stmcture of A 3 C6o- 
subspace that contains the degrees of freedom that we are interested in (in 
our case the 'ti u -subspace') and the rest of the Hilbert space: 7i = 7io®7i\. 
We can then write the Hamiltonian of the system as 



rr _ ( H 00 \ / V 01 \ m 

H -{ H n ) + [v w J' W 

where Ha is the projection of the Hamiltonian onto subspace TCi, while 
Vij = Hij contain the hybridization matrix elements between the two sub- 
spaces. Writing Green's function G = (E — H)~ l in the same way, we can 
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Figure 2. Band structure (ti„ band) of solid Ceo (fee) (a) as calculated ab initio using 
the local density approximation [1] and (b) using a tight-binding Hamiltonian with only 
tin orbitals [4]. 



calculate the projection of G onto TLq [3]: 

Goo = (E- [Hoo + V 01 (E - tfii)-V 10 ] 



(2) 



H cS (E) 



We see that the physics of the full system is described by an effective 
Hamiltonian H e g(E) that operates on the subspace TCq only. This drastic 
simplification comes, however, at a price: the effective Hamiltonian is energy 
dependent. In practice one approximates it with an energy-independent 
Hamiltonian H e s(Eo). This works well if we are only interested in energies 
close to Eq. In solid Cgo we have the fortunate situation that the bands 
retain the character of the molecular orbitals, since the hybridization matrix 
elements are small compared to the energy separations of the orbitals. In 
fact we can neglect the other bands altogether and get the hopping matrix 
elements t in j n r between the t\ u orbitals n and n' on molecules i and j 
directly from a tight-binding parameterization [4, 5]. Figure 2 shows the 
comparison of the ab initio t\ u band structure with the band structure 
obtained from the tight-binding Hamiltonian with only t± u orbitals. 

To get a realistic description of the electrons in the t\ u band we have to 
include the correlation effects which come from the Coulomb repulsion of 
electrons in t± u orbitals on the same molecule. The resulting Hamiltonian 
which describes the interplay of the hopping of electrons and their Coulomb 
repulsion has the form 



H 
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(3) 
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The on-site Coulomb interaction U can be calculated within density func- 
tional theory [6]. It is given by the increase in the energy of the t\ u level 
per electron that is added to one molecule of the system. It is important 
to avoid double counting in the calculation of U. While the relaxation of 
the occupied orbitals and the polarization of neighboring molecules has to 
be included in the calculation, excitations within the t\ u band must be 
excluded, since they are contained explicitly in the Hamiltonian (3). The 
results are consistent with experimental estimates [7, 8]: U ~ 1.2 — 1.4 eV. 
For comparison, the width of the t\ u band is in the range W « 0.5—0.85 eV. 

3. Quantum Monte Carlo 

We now turn to the question of how to calculate the ground state of a 
lattice Hamiltonian like (3). To simplify the notation most examples in the 
present section are for the simple Hubbard model (only one orbital per site, 
next neighbor hopping matrix elements Uj = —t) on a 2 dimensional square 
lattice: 

H = -tJ2 4 c i + U Y, n m- ( 4 ) 

The band width for this model is W = 8t. 

We first introduce the Gutzwiller Ansatz as a suitable trial function 
for the above Hamiltonian. Expectation values for the Gutzwiller function 
can be calculated using variational Monte Carlo (VMC). Then we describe 
the fixed-node diffusion Monte Carlo (FN-DMC) method that allows us to 
calculate more accurate variational estimates of the ground state energy 
(see the lecture notes by G. Bachelet for a more complete discussion of FN- 
DMC). The main emphasis of our discussion will be on the optimization of 
the trial function both in variational and fixed-node diffusion Monte Carlo. 

3.1. VARIATIONAL MONTE CARLO 

A good trial function for the Hubbard model has to balance the opposing 
tendencies of the hopping term and the interaction term: Without interac- 
tion (i.e. for U = 0) the ground state of the Hamiltonian (4) is the Slater 
determinant $ that maximizes the kinetic energy. Without hopping (t = 0) 
the interaction is minimized. Since only doubly occupied sites, i.e. sites with 
rijf = 1 and = 1, contribute to the Coulomb energy, the electrons are 
distributed as uniformly as possible over the lattice to minimize the number 
of double occupancies. A good compromise between these two extremes is 
to start from the non-interacting wavefunction $ but reduce the weight of 
configurations R with large double occupancies D(R). This leads (up to 
normalization) to the Gutzwiller wavefunction [9]: 



^ T {R) = g D{R) *(R), 



(5) 
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Figure 3. Weight of configurations with given number D of double occupancies for 
Gutzwiller wavefunctions ^t(R) = g D( - R ^ &(R). Reducing the Gutzwiller factor g sup- 
presses configurations with high Coulomb energy Ecoui{R) = U D(R) at the expense 
of increasing the kinetic energy. The results shown here are for a Hubbard model with 
16 x 16 sites and 101 + 101 electrons. 



with g £ (0, 1] the Gutzwiller parameter. Figure 3 shows how decreasing 
the Gutzwiller factor suppresses the configurations with a large number of 
double occupancies. 

To calculate the energy expectation value for the Gutzwiller wavefunc- 
tion we have to perform a sum over all configurations R: 



1 <*t|*t> E R n(R) 



where we have introduced the local energy for a configuration R 

Since the number of configurations R grows exponentially with system-size, 
the summation in (6) can be performed only for very small systems. For 
larger problems we use variational Monte Carlo [10]. The idea is to perform 
a random walk in the space of configurations, with transition probabilities 
p(R — > R') chosen such that the configurations Rvmc hi the random walk 
have the probability distribution function 9j,[R). Then 



EvMC= — ss« =Et - (8) 



The transition probabilities can be determined from detailed balance 

^l(R) p{R -» R!) = ^(R!) p{R! R) (9) 



which gives p(R -> R') = l/N min[l, V%,(R')/V%,(R)], with N being the 
maximum number of possible transitions. It is sufficient to consider only 
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transitions between configurations that are connected by the Hamiltonian, 
i.e. transitions in which one electron hops to a neighboring site. The stan- 
dard prescription is then to propose a transition R — ► R' with probability 
l/N and accept it with probability min[l, ^(R')/^(R)]. This works well 
for U not too large. For strongly correlated systems, however, the random 
walk will stay for long times in configurations with a small number of dou- 
ble occupancies D(R), since most of the proposed moves will increase D 
and hence be rejected with probability ~ 1 — g D<yR )~ D ( R \ 

Fortunately there is a way to integrate-out the time the walk stays in a 
given configuration. To see how, we first observe that for the local energy (7) 
the ratio of the wavefunctions for all transitions induced by the Hamiltonian 
have to be calculated. This in turn means that we also know all transition 
probabilities p(R — ► R'). We can therefore eliminate any rejection (i. e. 
accept with probability one) by proposing moves with probabilities 

p(R ^ R) -j: Rl p(R^R')-i-P^(Ry (10) 

Checking detailed balance (9) we find that now we are sampling configura- 
tions Rvmc from the probability distribution function ^^(i?) (l—p stay (R)). 
To compensate for this we assign a weight w(R) = 1/(1 — Pst& y (R)) to each 
configuration R. The energy expectation value is then given by 

Et w ^Rvmc w ^vmc) E 1oc {R VM c) 

^Rvmc W (RVMC) 

The above method is quite efficient since it ensures that in every Monte 
Carlo step a new configuration is created. Instead of staying in a configu- 
ration where ^>t is large, this configuration is weighted with the expecta- 
tion value of the number of steps the simple Metropolis algorithm would 
stay there. This is particularly convenient for simulations of systems with 
strong correlations: Instead of having to do longer and longer runs as U is 
increased, the above method produces, for a fixed number of Monte Carlo 
steps, results with comparable error estimates. 

Correlated sampling 

We now turn to the problem of optimizing the trial function *&t- A crite- 
rion for a good trial function is e.g. a low variational energy. To find the 
wavefunction that minimizes the variational energy we could do indepen- 
dent VMC calculations for a set of different trial functions. It is, however, 
difficult to compare the energies from these calculations since each VMC 
result comes with its own statistical errors. This problem can be avoided 
with correlated sampling [11]. The idea is to use the same random walk for 
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Figure 4- Correlated sampling for the 
Gutzwiller parameter g. The calculations 
are for a Hubbard model with 8x8 sites, 
28 + 28 electrons, and U = At. The full 
curve shows the result starting from a cal- 
culation with g = 1. The predicted mini- 
mum <? m j n is indicated by the dotted line. 
A dashed line gives the correlated sam- 
pling curve obtained from a calculation us- 
ing g m i n in the trial function. Both find the 
same minimum. E(g) becomes unreliable 
for very small g due to reweighting factors 
much larger than unity. 



calculating the expectation value of all the different trial functions. This 
reduces the relative errors and hence makes it easier to find the minimum. 

Let us assume we have generated a random walk Rvmc using as 
a trial function. Using the same random walk, we can then estimate the 
energy expectation value (8) for a different trial function by introducing 
the reweig hting factors ^%{R)/^{R): 

& ^ E RvMC E loc (R)mR)/n(R) (u) 

(Since the random walk Rvmc has only a finite number of configurations, 
this will only work well as long as the reweighting factors do not deviate 
too much from unity. Otherwise a few configurations with large reweight- 
ing factors will dominate. See Fig. 4.) We notice that (also in E\ oc ) the 
new trial function *$>t appears only in ratios with the old trial function. For 
Gutzwiller functions (5) that differ only in the Gutzwiller factor this means 
that the Slater determinants cancel, leaving only powers (g/g) D ^- Since 
D(R) is integer we can then rearrange the sums in (12) into polynomials 
in g/g. To find the optimal Gutzwiller parameter we then pick a reason- 
able g, perform a VMC run for $>T(g) during which we also estimate the 
coefficients for these polynomials. We can then calculate E(g) by simply 
evaluating the ratio of the polynomials. Since there are typically only of the 
order of some ten non-vanishing coefficients (cf. the distribution of weights 
in Fig. 3), this method is very efficient. Figure 4 shows how the method 
performs in practice. The idea of rewriting the sum over configurations into 
a polynomial can be easily generalized to trial functions with more corre- 
lation factors of the type r c( - R ^ , as long as the correlation function c(R) is 
integer-valued on the space of configurations. 
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Figure 5. Dependence of variational 
(VMC) and fixed-node diffusion Monte 
Carlo (FN-DMC) on the trial function. Uo 
is the Hubbard interaction that was used 
for the Slater determinant in the Gutzwiller 
wavefunction # T (i?) = g D(R) <f>(U ). The 
Gutzwiller parameter has always been op- 
timized. The results shown here are the 
energies (relative to the atomic limit) for 
a Hamiltonian that describes K3C60 (32 
sites), with U being varied from 1.25 (low- 
est curve) to 2.00 eV (highest curve). 

Character of the Slater determinant 
So far we have always constructed the Gutzwiller wavefunction from the 
ground state $ of the non-interacting Hamiltonian (U = 0). Alternatively 
we could use the Slater determinant &(U) from solving the interacting prob- 
lem in the Hartree-Fock approximation. We can even interpolate between 
these two extremes by doing a Hartree-Fock calculation with a fictitious 
Hubbard interaction Uq to obtain the Slater determinant $(Uq). This intro- 
duces an additional variational parameter in the Gutzwiller wavefunction. 
Increasing Uq will change the character of the trial function from paramag- 
netic to antiferromagnetic. This transition is also reflected in the variational 
energies, as is shown in Figure 5. Clearly, for small U the paramagnetic state 
is favorable, while for large U the antiferromagnetic state gives a lower vari- 
ational energy. We notice that for all values of U the optimal Uo is much 
smaller than U. 

3.2. FIXED-NODE DIFFUSION MONTE CARLO 

Diffusion Monte Carlo [12] allows us, in principle, to sample the true ground 
state of a Hamiltonian. The basic idea is to use a projection operator that 
has the lowest eigenstate as a fixed point. For a lattice problem where the 
spectrum is bounded E n G [Eo, E max ], the projection is given by 

|gr(n+i)) = [i- T (H-Eo)] |* (n) ) ; |* (0) ) = |¥ T >- (13) 

If r < 2/(E max — Eo) and \ ^t) has a non- vanishing overlap with the ground 
state, the above iteration converges to |^o)- There is no time-step error 
involved. Because of the prohibitively large dimension of the many-body 
Hilbert space, the matrix- vector product in (13) cannot be done exactly. 
Instead, we rewrite the equation in configuration space 

Y \R')(R'\y( n+ V) = Y \R') (R'\l - t{H - E )\R)(R\V {n) ) (14) 

=:F(R',R) 
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and perform the propagation in a stochastic sense: \E'( n ) is represented by 
an ensemble of configurations R with weights w(R). The transition ma- 
trix element F(R',R) is rewritten as a transition probability p(R — ► R') 
times a normalization factor m(R',R). The iteration (14) is then stochas- 
tically performed as follows: For each R we pick a new configuration R' 
with probability p(R — ► i?') and multiply its weight by m(R',R). Then 
the new ensemble of configurations i?' with their respective weights repre- 
sents \|/( n+1 ). Importance sampling decisively improves the efficiency of this 
process by replacing F(R',R) with G(R',R) = {<&t\R') F(R',R)/{R\V T ), 
so that transitions from configurations where the trial function is small to 
configurations with large trial function are enhanced: 

J2\R')(^ T \R')(R'\^ {n+1) ) = \R') G{R' , R) (V T \R} (R\V {n) }- (15) 

Now the ensemble of configurations represents the product \E , T 1 I/( n ). After 
a large number n of iterations the ground state energy is then given by the 
mixed estimator 

Eo - (^R) e^T' (16) 

As long as the evolution operator has only non-negative matrix elements 
G(R',R), all weights w(R) will be positive. If, however, G has negative 
matrix elements there will be both configurations with positive and neg- 
ative weight. Their contributions to the estimator (16) tend to cancel so 
that eventually the statistical error dominates, rendering the simulation 
useless. This is the infamous sign problem. A straightforward way to get 
rid of the sign problem is to remove the offending matrix elements from the 
Hamiltonian, thus defining a new Hamiltonian H e g by 

<JM> = { W V) ifG( *f <0 < 17 > 

For each off-diagonal element (R'\H\R) that has been removed, a term is 
added to the diagonal: 

(R\H eS \R) = (R\H\R) + J2 Vt{R')(R'\H\R)/Vt{R). (18) 

R' 

This is the fixed-node approximation for lattice Hamiltonians introduced 
in Ref. [13]. i? e ff is by construction free of the sign problem and variational, 
i.e. Ef > E Q . The equality holds if ^ t (R')/^t(R) = ^o(R')/^o(R) for 
all R, R' with G(R',R) < 0. 
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Figure 6. Correlated sampling of 
the Gutzwiller parameter g in the 
trial function to optimize the ef- 
fective Hamiltonian in fixed-node 
diffusion Monte Carlo. The results 
shown are for a Hubbard model 
with 4x4 sites, 7 + 7 electrons, 
and U = 4i. The error bars are the 
FN-DMC energies for different val- 
ues of g, the lines through the er- 
ror bars are the corresponding cor- 
related sampling curves. 



Fixed-node diffusion Monte Carlo for a lattice Hamiltonian thus means 
that we choose a trial function from which we construct an effective Hamil- 
tonian and determine its ground state by diffusion Monte Carlo. Because 
of the variational property, we want to pick the such that Eq S is min- 
imized, i.e. we want to optimize the trial function, or, equivalently, the 
effective Hamiltonian. Also here we can use the concept of correlated sam- 
pling. For optimizing the Gutzwiller parameter g we can even exploit the 
idea of rewriting the correlated sampling sums into polynomials mg/g, that 
we already have introduced in VMC. There is, however, a problem arising 
from the fact that the weight of a given configuration in iteration n is 
given by the product w(R^) = Yli=i m(R^\ R^ l ~^). Each individual nor- 
malization factor m(R',R) can be written as a finite polynomial, but the 
order of the polynomial for w(R^) increases steadily with the number of 
iterations. It is therefore not practical to try to calculate the ever increas- 
ing number of coefficients for the correlated sampling function E^ n \g). But 
since we still can easily calculate the coefficients for the m(R',R), we may 
use them to evaluate E^ (g) in each iteration on a set of predefined values 
gi of the Gutzwiller parameter. Figure 6 shows an example. It is interesting 
to note that the Gutzwiller factor that minimizes Eyuc is usually not the 
optimum Gutzwiller factor for fixed- node DMC. 

As in VMC we can also vary the trial function by changing the character 
of the Slater determiant &(Uq). We again find that the change from a 
paramagnetic to an antiferromagnetic trial function is reflected in the fixed- 
node energies (see Fig. 5), the paramagnetic state being favored for small 
U, while the antiferromagnetic state gives the lower energy for large U. 

We want to use Monte Carlo methods to detect a Mott transition in the 
doped Fullerides. For this we anticipate that we need an accuracy of better 
than 0.025 eV . To get a feeling for the accuracy of variational and fixed- 
node diffusion Monte Carlo, using Gutzwiller trial functions, we compare 
the results of the QMC calculations with exact results. Since exact diago- 
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Table 1. Total energy (in eV) for a cluster of four Ceo molecules with 6 + 6 electrons in 
the ti u band (hopping parameters for K 3 C6o)- The results of variational and diffusion 
Monte Carlo are compared to the exact energy. 
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0.038 
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23.2732 


23.2757(10) 


0.003 


23.3125(6) 
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nalizations can only be done for small systems we consider a small cluster 
of 4 molecules. The results for different values of the Hubbard interaction 
U are shown in Table 1. We find that the FN-DMC error is about an order 
of magnitude smaller than the error in VMC. The typical FN-DMC error 
for our lattice model is typically some meV, which should be sufficient for 
the application at hand. 

4. Mott transition in doped Fullerides 

We now apply the quantum Monte Carlo methods described above to the 
Hamiltonian (3). Our aim is to understand the Mott transition in the 
integer-doped Fullerides A ra C6o- Here A stands for an alkali metal like K, 
Rb, or Cs. The criterion for the metal-insulator transition is the opening 
of the gap 

Eg = E(N+l)-2E(N) + E(N -1). (19) 

Density functional calculations predict that the doped Fullerides A n C6o 
with n = 1 . . . 5 are metals [1]. Only AeCeo is an insulator with a completely 
filled t\ u band. On the other hand, the strong Coulomb repulsion between 
two electrons on the same Cqo molecule, which is much larger than the 
width of the t\ u band, suggests that all integer-doped Fullerides should be 
Mott insulators. It has therefore been suggested that experimental samples 
of, say, the superconductor K3C60 are metallic only because they are non- 
stoichiometric, i.e. that they actually are K^_sCqo [8]. 

In a first step we investigate what consequences the degeneracy of the 
ii u -band has for the Mott transition in K^Ceo- The analysis is motivated 
by the following simple argument [14, 15]. In the limit of very large U 
we can estimate the energies needed to calculate the gap (19). For half 
filling, all molecules will have 3 electrons in the t\ u orbital (Fig. 7, top). 
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Hopping is strongly suppressed since it would increase the energy by U. 
Therefore, to leading order in t 2 /U, there will be no kinetic contribution to 



the total energy E(N). In contrast, the systems with 
N ± 1 electrons have an extra electron/hole that can 
hop without additional cost in Coulomb energy To 
estimate the kinetic energy we calculate the matrix 
element for the hopping of the extra charge against 
an antiferromagnetic background. Denoting the ini- 
tial state with extra charge on molecule i by |1), we 
find that the second moment (l|ff 2 |l) is given by the 
number of different possibilities for a next-neighbor 
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Figure 7. 

hop times the single electron hopping matrix element ar g ument - 
t squared. By inserting £)■ \ where \j) denotes the state with the ex- 
tra charge hopped from site i to site j, we find (l\H\j) = y/3t, since, with 
an antiferromagnetic background and degeneracy 3, there are 3 different 
ways an extra charge can hop to a neighboring molecule (Fig. 7, bottom). 
Thus, due to the 3- fold degeneracy, the hopping matrix element is enhanced 
by a factor y/3 compared to the single electron hopping matrix element t. 
For a single electron system the kinetic energy is of the order of half the 
band width W/2. The enhancement of the hopping matrix element in the 
many-body case suggests then that the kinetic energy for the extra charge 
is correspondingly enhanced. Inserting the energies into (19) we find that 
for the 3-fold degenerate system our simple argument predicts a gap 



Eg = U-V3W, 



(20) 



instead of E g = U — W in the non-degenerate case. Extrapolating to in- 
termediate U, it appears that the degeneracy shifts the Mott transition 
towards larger U. 

The above argument is, of course, not rigorous. First, it is not clear 
whether the result for E g that was obtained in the limit of large U can be 
extrapolated to intermediate U, where the Mott transition actually takes 
place. Also the analogy of the hopping in the many-body case with the hop- 
ping of a single electron is not rigorous, since the hopping of an extra charge 
against an antiferromagnetic background creates a string of flipped spins 
[16]. Nevertheless the argument suggests that orbital degeneracy might play 
an important role for the Mott transition. 

To test this proposition, we have performed quantum Monte Carlo cal- 
culations for the model Hamiltonian (3) with hopping matrix elements ap- 
propriate for K3C60 [14]. The Coulomb interaction U has been varied from 
U = . . . 1.75 eV to study the opening of the gap. Since the Monte Carlo 
calculations are for finite systems, we have to extrapolate to infinite sys- 
tem size. To improve the extrapolation we correct for finite-size effects: 
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Figure 8. Finite-size corrected gap 
E g = E g - U/M - E g {U = 0) for 
increasing Coulomb interaction U as 
a function of 1/M, where M is the 
number of molecules. The calculations 
are for a Hubbard model with hop- 
ping matrix elements appropriate for 
K^Ceo- The band width varies be- 
tween W = 0.58 eV for M = 4 and 
W — 0.63 eV in the infinite-size limit. 



First, there could be a gap E g (U = 0) already in the spectrum of the 
non-interacting system. Further, even for a metallic system of M molecules 
there will be a finite-size contribution of U/M to the gap. It comes from 
the electrostatic energy of the extra charge, uniformly distributed over all 
sites. Both corrections vanish in the limit M — > oo, as they should. The 
finite-size corrected gap E g = E g — U/M — E g (U = 0) for systems with 
M = 4, 8, 16, 32, and 64 molecules is shown in Figure 8. We find that the 
gap opens for U between 1.50 eV and 1.75 eV. Since for the real system 
U =1.2... 1.4 eV, K3C60 is thus close to a Mott transition, but still on the 
metallic side - even though U is considerably larger than the band width 
W. This is in contrast to simpler theories that neglect orbital degeneracy. 

Doping dependence 

The degeneracy argument described 
above for K3C60 can be generalized to 
integer fillings. Away from half filling 
the enhancement of the hopping ma- 
trix elements for an extra electron is 
different from that for an extra hole. 
The effective enhancement for different 
fillings are given in the adjacent table. 

We find that the enhancement decreases as we move away from half 
filling. Therefore we expect that away from half filling, correlations be- 
come more important, putting the system closer to the Mott transition, or 
maybe even pushing it across the transition, making it an insulator. We 
have analyzed the doping dependence of the Mott transition for the same 
Hamiltonian as used for K3C60, changing the filling of the t\ u band from 
n = 1 to 5 [17]. This model describes the Fm3m-Fullerides A n C6o with fee 
lattice and orientational disorder [18]. The critical Coulomb interaction U c , 
at which the transition from a metal (for U < U c ) to an insulator (U > U c ) 
takes place, is shown in Figure 9 for the different integer fillings. As ex- 
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Figure 9. Doping dependence of the 
Mott transition. The error bars indi- 
cate the estimate of the critical ratio 
Uc/W for different integer fillings of 
the tiu band. The calculations are for 
doped Fm3m Fullerides with fee lat- 
tice structure and orientational disor- 
der. The shaded region shows the range 
of U /W in which the doped Fullerides 
are falling. 




pected from the degeneracy argument, U c decreases away from n = 3. We 
note, however, that U c is asymmetric around half filling. This asymmetry 
is not present in the simple degeneracy argument, where we implicitly as- 
sumed that the lattice is bipartite. In such a situation we have electron-hole 
symmetry, which implies symmetry around half-filling. For frustrated lat- 
tices like the fee lattice electron-hole symmetry is broken, leading to the 
asymmetry in U c that is seen in Fig. 9. 

Lattice dependence 

To understand the effect of frustration in terms of the hopping argu- 
ments that we have made so far, we have to consider more 
than just one next-neighbor hop. The simplest system where 
we encounter frustration is a triangle with hopping between 
neighboring sites. In the single electron case we can form a # - # 
bonding state with energy E m i n = 2 1, but because of frustra- Figure 10. : 
tion we cannot form an antibonding state. Instead the maxi- Tnan gl e 
mum eigenenergy is E max = t. Hence frustration leads to an asymmetric 
'band' of width W = 3t. 

In the many-body case the situation is different. Like in the degeneracy 
argument we look at the hopping of an extra electron against a (frustrated) 
antiferromagnetic background in the large-?/ limit. For simplicity we assume 
a non-degenerate system, i.e. there is one electron per site on the triangle, 
plus the extra electron. In this case we have to move the extra charge 
twice around the triangle to come back to the many-body state we started 
from. Thus in the large-£7 limit the many-body problem is an eigenvalue 
problem of a 6 x 6 matrix with extreme eigenvalues ±2 1. In the degeneracy 
argument we have assumed that the kinetic energy of the extra charge is 
given by W/2. On the triangle, we find, however, that the hopping energy is 
larger than that by a factor 4/3. This suggests that for frustrated systems 
the single electron band width W in (20) should be multiplied by a prefactor 
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larger than one. We therefore expect that frustration alone, even without 
degeneracy, shifts the Mott transition to larger U. 

To analyze the effect of frustration on the Mott transition we have de- 
termined the critical U for a hypothetical doped Fulleride A4C60 with body 
centered tetragonal (bet) structure, a lattice without frustration, having the 
same band width (W = 0.6 eV) as the fcc-Fullerides, shown in Figure 9. For 
U = 1.3 eV , we find a gap E g ~ 0.6 eV for the Fulleride with bet structure, 
while the frustrated fee compound is still metallic E g = 0. This difference 
is entirely due to the lattice structure. Using realistic parameters for K4C60 
[5] that crystallizes in a bet structure we find a Mott insulator with gap 
Eg ~ 0.7 eV, which is in line with experimental findings: E g = 0.5 ± 0.1 eV 
[19]. 

Conclusion 

We have seen that, due to more efficient hopping, orbital degeneracy in- 
creases the critical U at which the Mott transition takes place. This puts 
the integer-doped Fullerides close to a Mott transition. Whether they are 
on the metallic or insulating side depends on the filling of the band and 
the lattice structure: Since the degeneracy enhancement works best for a 
half filled band, systems doped away from half-filling tend to be more insu- 
lating. The effect of frustration, on the other hand, is to make the system 
more metallic. 
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